# 导入所需的库
import pandas as pd  # 用于数据读取和处理，包括读取 Excel 数据文件
import matplotlib.pyplot as plt  # 用于创建可视化图形，支持 2D 和 3D 绘图
import seaborn as sns  # 高级绘图库，能够产生美观的统计图形，增强 matplotlib 的功能
from sklearn.mixture import GaussianMixture  # 使用高斯混合模型（GMM）进行聚类分析
from sklearn.preprocessing import StandardScaler  # 用于将数据标准化（均值为0，方差为1）
import numpy as np  # 数值计算库，提供多维数组和矩阵处理功能
from mpl_toolkits.mplot3d import Axes3D  # 用于创建 3D 图形的工具包

# 设置 Matplotlib 的全局配置
# 指定字体为 SimHei，解决中文显示问题（如标题、轴标签等）
plt.rcParams['font.sans-serif'] = ['SimHei']
# 防止坐标轴上的负号显示为方块，启用负号的正确显示
plt.rcParams['axes.unicode_minus'] = False
# 将默认字体大小设置为 14，确保绘图中字体大小的一致性
plt.rcParams['font.size'] = 14

# 读取 Excel 文件中的数据
file_path = '附件2.xlsx'  # 文件路径
# 从 Excel 中读取特定的工作表，这里选择的是“2023年统计的相关数据”表
data = pd.read_excel(file_path, sheet_name='2023年统计的相关数据')

# 由于“销售单价/(元/斤)”是一个范围值，需要将其转换为平均值
# 使用 lambda 函数将字符串格式的区间（如 '10-20'）转换为两个数字的平均值
data['平均销售单价'] = data['销售单价/(元/斤)'].apply(lambda x: sum(map(float, x.split('-'))) / 2)

# 提取与聚类分析相关的列：种植成本、亩产量、平均销售单价
# 将这些列用于后续的聚类分析
data_for_clustering = data[['种植成本/(元/亩)', '亩产量/斤', '平均销售单价']]

# 对数据进行标准化处理，使不同特征（如种植成本和亩产量）具有相同的尺度
# 标准化后的数据具有零均值和单位方差
scaler = StandardScaler()  # 初始化 StandardScaler 对象
data_scaled = scaler.fit_transform(data_for_clustering)  # 标准化数据，结果为 numpy 数组

# 使用高斯混合模型 (Gaussian Mixture Model, GMM) 进行聚类分析
# n_components=3 表示我们将数据聚成 3 类
gmm = GaussianMixture(n_components=3, random_state=0)  # 初始化 GMM 模型，假设 3 类
gmm.fit(data_scaled)  # 拟合标准化后的数据
clusters = gmm.predict(data_scaled)  # 预测每个数据点所属的聚类类别

# 将聚类结果添加到原数据中，作为新列“聚类类别”
# 每个数据点根据其特征被归入 0, 1, 或 2 类之一
data['聚类类别'] = clusters

# 创建三维绘图，展示不同聚类的农作物分布
fig = plt.figure(figsize=(10, 8))  # 设置绘图窗口的大小
ax = fig.add_subplot(111, projection='3d')  # 创建三维图形的坐标系

# 设置三维图的视角，elev 是从上往下的角度，azim 是从侧面旋转的角度
ax.view_init(elev=30, azim=30)

# 绘制三维散点图，X 轴表示种植成本，Y 轴表示亩产量，Z 轴表示平均销售单价
# 使用不同的颜色来表示不同的聚类类别
scatter = ax.scatter(
    data['种植成本/(元/亩)'], data['亩产量/斤'], data['平均销售单价'],  # X, Y, Z 轴数据
    c=data['聚类类别'], cmap='viridis', s=100, marker='o', edgecolor='k'  # 根据聚类结果着色
)

# 设置三维图形的标题和坐标轴标签
ax.set_title('农作物的三维聚类分析 (GMM)', fontsize=30, weight='bold', pad=25)  # 标题
ax.set_xlabel('种植成本/(元/亩)', fontsize=15)  # X 轴标签
ax.set_ylabel('亩产量/斤', fontsize=15)  # Y 轴标签
ax.set_zlabel('平均销售单价', fontsize=15)  # Z 轴标签

# 添加颜色条以指示不同的聚类类别
cbar = fig.colorbar(scatter, ax=ax, pad=0.1)  # 创建颜色条
cbar.set_label('聚类类别', fontsize=20)  # 设置颜色条标签

# 调整布局，使得子图之间没有重叠
plt.tight_layout()

# 显示绘制的三维图
plt.show()

# 函数：plot_ellipsoid
# 功能：根据给定的均值和协方差矩阵绘制置信椭球
# 椭球表示数据点的分布情况，置信椭球可以帮助可视化某个聚类中的数据分布范围
def plot_ellipsoid(mean, cov_matrix, ax):
    """
    根据均值和协方差矩阵绘制置信椭球。

    :param mean: 均值向量，表示椭球中心
    :param cov_matrix: 协方差矩阵，表示数据的分布形状
    :param ax: 三维绘图的轴对象，用于绘制椭球
    """
    # 通过奇异值分解 (SVD) 计算协方差矩阵的特征值和特征向量
    # U 是旋转矩阵，s 是特征值，rotation 是旋转矩阵
    U, s, rotation = np.linalg.svd(cov_matrix)
    radii = np.sqrt(s)  # 特征值的平方根对应椭球的半径

    # 生成网格点，用于绘制椭球的表面
    u = np.linspace(0, 2 * np.pi, 100)  # 参数 u，定义椭球的角度范围（0 到 2π）
    v = np.linspace(0, np.pi, 100)  # 参数 v，定义椭球的纵向角度范围（0 到 π）
    x = radii[0] * np.outer(np.cos(u), np.sin(v))  # 椭球的 x 坐标
    y = radii[1] * np.outer(np.sin(u), np.sin(v))  # 椭球的 y 坐标
    z = radii[2] * np.outer(np.ones_like(u), np.cos(v))  # 椭球的 z 坐标

    # 将网格点进行旋转和平移，使其符合数据的分布
    for i in range(len(x)):
        for j in range(len(x)):
            # 对每个网格点进行旋转，并将其移动到均值的位置
            [x[i, j], y[i, j], z[i, j]] = np.dot([x[i, j], y[i, j], z[i, j]], rotation) + mean

    # 绘制椭球的线框
    ax.plot_wireframe(x, y, z, color='r', alpha=0.1)  # 设置颜色和透明度

# 创建三维图形并绘制带有置信椭球的聚类图
fig = plt.figure(figsize=(10, 8))  # 创建新的绘图窗口
ax = fig.add_subplot(111, projection='3d')  # 创建三维图形的坐标系

# 绘制聚类点，X 轴表示种植成本，Y 轴表示亩产量，Z 轴表示平均销售单价
scatter = ax.scatter(
    data['种植成本/(元/亩)'], data['亩产量/斤'], data['平均销售单价'],  # X, Y, Z 轴数据
    c=data['聚类类别'], cmap='inferno', s=100, marker='o', edgecolor='k'  # 根据聚类结果着色
)

# 设置三维图形的标题和坐标轴标签
ax.set_title('农作物的置信椭球 (GMM)', fontsize=30, weight='bold', pad=20)  # 标题
ax.set_ylabel('亩产量/斤', fontsize=15)  # Y 轴标签
ax.set_xlabel('种植成本/(元/亩)', fontsize=15)  # X 轴标签
ax.set_zlabel('平均销售单价', fontsize=15)  # Z 轴标签

# 添加颜色条，表示不同的聚类类别
cbar = fig.colorbar(scatter, ax=ax, pad=0.1)
cbar.set_label('聚类类别', fontsize=20)

# 为每个聚类添加置信椭球，帮助可视化各类的分布形状
for i in range(3):  # 遍历 3 个聚类
    # 获取当前聚类的数据点
    cluster_points = data_for_clustering[clusters == i]

    # 如果当前聚类点的数量大于 2，则绘制置信椭球
    if len(cluster_points) > 2:
        # 计算当前聚类的均值和协方差矩阵
        mean = cluster_points.mean(axis=0)
        cov_matrix = np.cov(cluster_points.T)

        # 绘制该聚类的置信椭球
        plot_ellipsoid(mean, cov_matrix, ax)

# 调整图形布局，确保没有重叠
plt.tight_layout()

# 显示图形
plt.show()
